# =============================================================================
# Statistical Analysis Code for:
# "Chronic adaptive versus conventional DBS response patterns in Parkinson's
#  disease: A pilot randomized crossover trial"
#
# Main Analysis: Standard Crossover ANCOVA with Mixed Effects Models
# Author: Jun Tanimura, MD, MSc.
# Created: 2025-08-07
# Last Updated: 2025-11-26
# Version: v3.1
#
# CHANGES FROM v3.0 FOR PIEER REVIEW REVISION:
# - Added 95% confidence intervals for Cohen's d effect sizes
# - CI calculated using t-distribution (consistent with treatment effect CI)
# =============================================================================

library(tidyverse)
library(lme4)
library(lmerTest)
library(car)
library(broom)
library(emmeans)

# Read the data
df <- read_csv("data.csv")

# Define evaluation items in order
eval_items_ordered <- c("ADL", "UPDRS_part_1", "UPDRS_part_2", "UPDRS_part_3", "UPDRS_part_4", 
                        "UPDRS_total", "Mean_corr_ON_time", "Mean_corr_OFF_time", 
                        "Mean_corr_Dys_sev_time", "Mean_corr_Dys_weak_time", "MMSE", "PDQ-39")

# Function to adjust Cohen's d sign for directional interpretation
adjust_cohens_d_for_direction <- function(metric_name, raw_effect_size) {
  # Define metrics where lower values are better
  lower_the_better_metrics <- c(
    "UPDRS_part_1", "UPDRS_part_2", "UPDRS_part_3", "UPDRS_part_4", "UPDRS_total",
    "PDQ-39",
    "Mean_corr_Dys_sev_time", "Mean_corr_Dys_weak_time",
    "Mean_corr_OFF_time"
  )
  
  # For lower-the-better metrics, reverse the sign (negative treatment effect = improvement -> positive Cohen's d)
  if (metric_name %in% lower_the_better_metrics) {
    return(-raw_effect_size)
  } else {
    # Higher the better (ON_time, MMSE) - keep as is
    return(raw_effect_size)
  }
}

# Function to test random intercept significance using LRT
test_random_intercept_lrt <- function(complete_data) {
  cat("Testing random intercept significance with LRT...\n")
  
  # Model 0: No random intercept (fixed effects only)
  m0_ml <- tryCatch({
    lm(Score_transformed ~ Treatment + Period + Sequence + Baseline_z, 
       data = complete_data)
  }, error = function(e) {
    cat(sprintf("Fixed effects model failed: %s\n", e$message))
    return(NULL)
  })
  
  # Model 1: Random intercept included (same fixed effects + (1|PatientID))
  m1_ml <- tryCatch({
    lmer(Score_transformed ~ Treatment + Period + Sequence + Baseline_z + (1|PatientID), 
         data = complete_data, REML = FALSE, control = lmerControl(optimizer = "bobyqa"))
  }, error = function(e) {
    cat(sprintf("Random intercept model failed: %s\n", e$message))
    return(NULL)
  })
  
  if (is.null(m0_ml) || is.null(m1_ml)) {
    cat("Random intercept test models failed to fit\n")
    return(list(
      ri_lrt_p = NA_real_,
      ri_lrt_stat = NA_real_,
      ri_var_subject = NA_real_,
      ri_var_residual = NA_real_,
      ri_icc = NA_real_
    ))
  }
  
  # Calculate LRT statistic
  ll0 <- logLik(m0_ml)
  ll1 <- logLik(m1_ml)
  lr_stat <- as.numeric(2 * (ll1 - ll0))
  
  # Boundary hypothesis p-value: 0.5 * (1 - pchisq(LR, df=1))
  # Uses mixture distribution since random effect variance = 0 is on the boundary
  ri_p <- 0.5 * (1 - pchisq(lr_stat, df = 1))
  
  # Variance components and ICC (from REML model)
  m1_reml <- tryCatch({
    lmer(Score_transformed ~ Treatment + Period + Sequence + Baseline_z + (1|PatientID), 
         data = complete_data, REML = TRUE, control = lmerControl(optimizer = "bobyqa"))
  }, error = function(e) {
    cat("REML model for variance components failed\n")
    return(NULL)
  })
  
  if (!is.null(m1_reml)) {
    vc <- as.data.frame(VarCorr(m1_reml))
    var_subj <- vc$vcov[vc$grp == "PatientID"]
    var_res <- sigma(m1_reml)^2
    icc <- var_subj / (var_subj + var_res)
  } else {
    var_subj <- NA_real_
    var_res <- NA_real_
    icc <- NA_real_
  }
  
  cat(sprintf("Random intercept LRT: LR = %.4f, p = %.4f\n", lr_stat, ri_p))
  cat(sprintf("Variance components: Subject = %.4f, Residual = %.4f, ICC = %.4f\n", 
              var_subj, var_res, icc))
  
  return(list(
    ri_lrt_p = ri_p,
    ri_lrt_stat = lr_stat,
    ri_var_subject = var_subj,
    ri_var_residual = var_res,
    ri_icc = icc
  ))
}

# Function to calculate time difference CI using Delta method
calculate_time_difference_delta <- function(model, metric_name, model_type = "mixed") {
  if (!str_detect(metric_name, "time")) {
    return(list(
      time_difference = NA,
      time_diff_ci_lower = NA,
      time_diff_ci_upper = NA,
      aDBS_time_adjusted = NA,
      cDBS_time_adjusted = NA,
      intuitive_ratio = NA,
      percent_reduction = NA,
      ratio_ci_lower = NA,
      ratio_ci_upper = NA
    ))
  }
  
  tryCatch({
    # Get marginal estimates and variance-covariance matrix from emmeans
    emm <- emmeans(model, "Treatment")
    emm_summary <- summary(emm)
    vcov_emm <- vcov(emm)
    
    # Check treatment order and get data
    treatment_levels <- emm_summary$Treatment
    aDBS_idx <- which(treatment_levels == "aDBS")
    cDBS_idx <- which(treatment_levels == "cDBS")
    
    if (length(aDBS_idx) == 0 || length(cDBS_idx) == 0) {
      warning("Treatment levels not found in emmeans results")
      return(list(
        time_difference = NA, time_diff_ci_lower = NA, time_diff_ci_upper = NA,
        aDBS_time_adjusted = NA, cDBS_time_adjusted = NA,
        intuitive_ratio = NA, percent_reduction = NA,
        ratio_ci_lower = NA, ratio_ci_upper = NA
      ))
    }
    
    # Estimates on log1p scale
    log1p_aDBS <- emm_summary$emmean[aDBS_idx]
    log1p_cDBS <- emm_summary$emmean[cDBS_idx]
    
    # Estimates on log1p scale
    aDBS_time <- exp(log1p_aDBS) - 1
    cDBS_time <- exp(log1p_cDBS) - 1
    
    time_difference <- aDBS_time - cDBS_time
    
    # Back-transform to time scale: exp(log1p_value) - 1
    grad_aDBS <- exp(log1p_aDBS)
    grad_cDBS <- exp(log1p_cDBS)
    
    # Gradient vector for time difference (aDBS - cDBS)
    gradient <- numeric(length(treatment_levels))
    gradient[aDBS_idx] <- grad_aDBS   # aDBS ã®ä¿‚æ•°
    gradient[cDBS_idx] <- -grad_cDBS  # cDBS ã®ä¿‚æ•°ï¼ˆå·®ã‚’å–ã‚‹ã®ã§è² ç¬¦å·ï¼‰
    
    var_time_diff <- as.numeric(t(gradient) %*% vcov_emm %*% gradient)
    se_time_diff <- sqrt(var_time_diff)
    
    # CI calculation - use Satterthwaite df correctly
    df_aDBS <- emm_summary$df[aDBS_idx]
    df_cDBS <- emm_summary$df[cDBS_idx]
    df_combined <- min(df_aDBS, df_cDBS)
    t_val <- qt(0.975, df_combined)
    
    time_diff_ci_lower <- time_difference - t_val * se_time_diff
    time_diff_ci_upper <- time_difference + t_val * se_time_diff
    
    ratio <- aDBS_time / cDBS_time
    percent_reduction <- (1 - ratio) * 100
    
  # Boundary hypothesis p-value: 0.5 * (1 - pchisq(LR, df=1))
    aDBS_ci_lower <- exp(emm_summary$lower.CL[aDBS_idx]) - 1
    aDBS_ci_upper <- exp(emm_summary$upper.CL[aDBS_idx]) - 1
    cDBS_ci_lower <- exp(emm_summary$lower.CL[cDBS_idx]) - 1
    cDBS_ci_upper <- exp(emm_summary$upper.CL[cDBS_idx]) - 1
    
    ratio_ci_lower <- aDBS_ci_lower / cDBS_ci_upper
    ratio_ci_upper <- aDBS_ci_upper / cDBS_ci_lower
    
    return(list(
      time_difference = time_difference,
      time_diff_ci_lower = time_diff_ci_lower,
      time_diff_ci_upper = time_diff_ci_upper,
      aDBS_time_adjusted = aDBS_time,
      cDBS_time_adjusted = cDBS_time,
      intuitive_ratio = ratio,
      percent_reduction = percent_reduction,
      ratio_ci_lower = ratio_ci_lower,
      ratio_ci_upper = ratio_ci_upper
    ))
    
  }, error = function(e) {
    warning(paste("Delta method calculation failed for", metric_name, ":", e$message))
    return(list(
      time_difference = NA, time_diff_ci_lower = NA, time_diff_ci_upper = NA,
      aDBS_time_adjusted = NA, cDBS_time_adjusted = NA,
      intuitive_ratio = NA, percent_reduction = NA,
      ratio_ci_lower = NA, ratio_ci_upper = NA
    ))
  })
}

# Transform data to long format for ANCOVA
ancova_data <- map_dfr(eval_items_ordered, function(item) {
  col_P <- paste0("P_", item)
  col_C <- paste0("C_", item)
  col_A <- paste0("A_", item)
  
  if(all(c(col_P, col_C, col_A) %in% names(df))) {
    # Create long format data
    patient_data <- df %>%
      select(PatientID, Group, all_of(c(col_P, col_C, col_A))) %>%
      pivot_longer(
        cols = c(all_of(col_C), all_of(col_A)),
        names_to = "Treatment_raw",
        values_to = "Score_post"
      ) %>%
      mutate(
        Baseline = !!sym(col_P),
        Treatment = case_when(
          str_detect(Treatment_raw, "^C_") ~ "cDBS",
          str_detect(Treatment_raw, "^A_") ~ "aDBS"
        ),
        Period = case_when(
          Group == 1 & Treatment == "cDBS" ~ "Phase1",
          Group == 1 & Treatment == "aDBS" ~ "Phase2", 
          Group == 2 & Treatment == "aDBS" ~ "Phase1",
          Group == 2 & Treatment == "cDBS" ~ "Phase2"
        ),
        Evaluation = item
      ) %>%
      mutate(
        # Apply log1p transformation for time variables
        is_time_var = str_detect(item, "time"),
        Score_transformed = ifelse(is_time_var, log1p(pmax(Score_post, 0)), Score_post),
        Baseline_transformed = ifelse(is_time_var, log1p(pmax(Baseline, 0)), Baseline)
      ) %>%
      select(PatientID, Group, Evaluation, Treatment, Period, 
             Baseline, Score_post, Baseline_transformed, Score_transformed)
  }
}) %>%
  filter(!is.na(Score_post)) %>%
  mutate(
    Treatment = factor(Treatment, levels = c("cDBS", "aDBS")),
    Period = factor(Period, levels = c("Phase1", "Phase2")),
    Evaluation = factor(Evaluation, levels = eval_items_ordered),
    PatientID = factor(PatientID)
  )

# Add z-score normalization for baseline values within each evaluation
ancova_data <- ancova_data %>%
  group_by(Evaluation) %>%
  mutate(
    Baseline_z = as.numeric(scale(Baseline_transformed))
  ) %>%
  ungroup()

print("=== Data transformation completed ===")
print(paste("Total observations:", nrow(ancova_data)))
print(paste("Unique evaluations:", length(unique(ancova_data$Evaluation))))

# Function to run standard crossover ANCOVA analysis (Updated with Random Intercept LRT AND Delta Method)
run_mixed_effects_analysis <- function(data, metric) {
  metric_data <- data %>% filter(Evaluation == metric)
  
  cat(sprintf("\n=== Standard Crossover ANCOVA Analysis for %s ===\n", metric))
  cat(sprintf("N observations: %d\n", nrow(metric_data)))
  
  # Check for missing values
  complete_data <- metric_data %>% 
    filter(!is.na(Score_transformed) & !is.na(Baseline_z))
  
  if(nrow(complete_data) < 10) {
    cat("Insufficient data for analysis\n")
    return(NULL)
  }
  
  # Add sequence variable (Group in our case)
  complete_data <- complete_data %>%
    mutate(Sequence = factor(Group))  # Group 1: cDBS->aDBS, Group 2: aDBS->cDBS
  
  # STEP 1: Test for differential carryover effect via LRT (minimal model: Treatment + Period only)
  cat("Testing differential carryover with minimal model (Treatment + Period + random effects only)...\n")
  
  m0_ml <- tryCatch({
    lmer(Score_transformed ~ Treatment + Period + (1|PatientID), 
         data = complete_data, REML = FALSE, control = lmerControl(optimizer = "bobyqa"))
  }, error = function(e) {
    cat(sprintf("Null model failed: %s\n", e$message))
    return(NULL)
  })
  
  m1_ml <- tryCatch({
    lmer(Score_transformed ~ Treatment * Period + (1|PatientID), 
         data = complete_data, REML = FALSE, control = lmerControl(optimizer = "bobyqa"))
  }, error = function(e) {
    cat(sprintf("Interaction model failed: %s\n", e$message))
    return(NULL)
  })
  
  if (is.null(m0_ml) || is.null(m1_ml)) {
    cat("Carryover test models failed to fit\n")
    return(NULL)
  }
  
  carryover_test <- tryCatch(anova(m0_ml, m1_ml), error = function(e) {
    cat(sprintf("ANOVA comparison failed: %s\n", e$message))
    return(NULL)
  })
  
  carryover_p <- if (!is.null(carryover_test)) carryover_test[["Pr(>Chisq)"]][2] else NA_real_
  
  # Clear reporting of carryover test results
  msg_p <- if (is.na(carryover_p)) "not computed" else sprintf("%.4f", carryover_p)
  cat(sprintf("Differential carryover test (LRT, minimal model): p = %s\n", msg_p))
  
  # STEP 2: Random Intercept LRT Test
  ri_test_results <- test_random_intercept_lrt(complete_data)
  
  # Always use full crossover analysis
  primary_analysis_type <- "Full_Crossover"
  
  # Fit standard crossover model without interaction (REML for primary estimation with baseline)
  primary_model <- tryCatch({
    lmer(Score_transformed ~ Treatment + Period + Sequence + Baseline_z + (1|PatientID), 
         data = complete_data, REML = TRUE, control = lmerControl(optimizer = "bobyqa"))
  }, error = function(e) {
    cat("Primary crossover model (REML) failed:", e$message, "\n")
    return(NULL)
  })
  
  analysis_data <- complete_data
  
  if(is.null(primary_model)) {
    cat("Primary analysis model failed\n")
    return(NULL)
  }
  
  # Extract results from primary model
  model_summary <- summary(primary_model)
  
  # Check if treatment effect exists
  coef_names <- rownames(model_summary$coefficients)
  treatment_coef_name <- "TreatmentaDBS"
  
  if(!treatment_coef_name %in% coef_names) {
    cat("Treatment coefficient not found in primary model\n")
    return(NULL)
  }
  
  treatment_effect <- model_summary$coefficients[treatment_coef_name, ]
  
  # Calculate effect size and other metrics
  residual_sd <- sigma(primary_model)
  # Calculate ICC for mixed models
  variance_components <- as.data.frame(VarCorr(primary_model))
  random_var <- variance_components$vcov[variance_components$grp == "PatientID"]
  residual_var <- variance_components$vcov[variance_components$grp == "Residual"]
  icc <- random_var / (random_var + residual_var)
  treatment_t <- treatment_effect[4]  # t value
  treatment_p <- treatment_effect[5]  # Satterthwaite p-value
  
  # Calculate Cohen's D and apply direction adjustment
  effect_size_raw <- treatment_effect[1] / residual_sd
  effect_size <- adjust_cohens_d_for_direction(metric, effect_size_raw)
  
  # Calculate 95% CI for Cohen's D (using accurate standard error formula)
  # Use Satterthwaite df
  cohens_d_df <- treatment_effect[3]
  
  # Accurate SE for Cohen's d in paired/crossover design
  # SE(d) = sqrt(1/n + d²/(2*df))
  n_patients <- length(unique(analysis_data$PatientID))
  se_cohens_d <- sqrt(1/n_patients + (effect_size_raw^2)/(2*cohens_d_df))
  t_critical <- qt(0.975, cohens_d_df)
  
  # Raw Cohen's d CI (before direction adjustment)
  ci_lower_d_raw <- effect_size_raw - t_critical * se_cohens_d
  ci_upper_d_raw <- effect_size_raw + t_critical * se_cohens_d
  
  # Cohen's d CI after direction adjustment
  ci_lower_d <- adjust_cohens_d_for_direction(metric, ci_lower_d_raw)
  ci_upper_d <- adjust_cohens_d_for_direction(metric, ci_upper_d_raw)
  
  # Check CI order (may be reversed by direction adjustment)
  if(ci_lower_d > ci_upper_d) {
    temp <- ci_lower_d
    ci_lower_d <- ci_upper_d
    ci_upper_d <- temp
  }
  
  ci_result <- tryCatch({
    emm_treatment <- emmeans(primary_model, "Treatment")
    emm_contrast <- contrast(emm_treatment, method = list("aDBS - cDBS" = c(-1, 1)))
    ci_result <- confint(emm_contrast)
    list(lower = ci_result$lower.CL[1], upper = ci_result$upper.CL[1])
  }, error = function(e) {
    cat("Confidence interval calculation failed\n")
    list(lower = NA, upper = NA)
  })
  
  # Calculate time difference using Delta method
  time_delta_results <- calculate_time_difference_delta(primary_model, metric, primary_analysis_type)
  
  # Print results
  cat("\n=== PRIMARY ANALYSIS RESULTS ===\n")
  cat(sprintf("Analysis type: %s\n", primary_analysis_type))
  cat(sprintf("Sample size: %d observations from %d patients\n", 
              nrow(analysis_data), length(unique(analysis_data$PatientID))))
  cat(sprintf("aDBS vs cDBS estimate: %.4f\n", treatment_effect[1]))
  cat(sprintf("Standard error: %.4f\n", treatment_effect[2]))
  cat(sprintf("t-value: %.4f\n", treatment_t))  # ä¿®æ­£: æ­£ã—ã„tå€¤ã‚’è¡¨ç¤º
  cat(sprintf("p-value: %.4f\n", treatment_p))
  cat(sprintf("Effect size (Cohen's d, aDBS-favorable direction): %.4f [95%% CI: %.4f, %.4f]\n", effect_size, ci_lower_d, ci_upper_d))
  if(!is.na(icc)) {
    cat(sprintf("ICC: %.4f\n", icc))
  }
  cat(sprintf("95%% CI: [%.4f, %.4f]\n", ci_result$lower, ci_result$upper))
  
  if (str_detect(metric, "time") && !is.na(time_delta_results$time_difference)) {
    cat("\n=== TIME DIFFERENCE RESULTS (Delta Method) ===\n")
    cat(sprintf("aDBS adjusted time: %.2f hours\n", time_delta_results$aDBS_time_adjusted))
    cat(sprintf("cDBS adjusted time: %.2f hours\n", time_delta_results$cDBS_time_adjusted))
    cat(sprintf("Time difference (aDBS - cDBS): %.2f hours\n", time_delta_results$time_difference))
    cat(sprintf("Time difference 95%% CI: [%.2f, %.2f] hours\n", 
                time_delta_results$time_diff_ci_lower, time_delta_results$time_diff_ci_upper))
    cat(sprintf("Intuitive ratio (aDBS/cDBS): %.4f\n", time_delta_results$intuitive_ratio))
    cat(sprintf("Percent change: %.1f%%\n", time_delta_results$percent_reduction))
  }
  
  # Return results (INCLUDING Random Intercept LRT results AND time delta)
  result_list <- list(
    metric = metric,
    analysis_type = primary_analysis_type,
    carryover_p = carryover_p,
    # Random Intercept LRT results (NEW)
    ri_lrt_p = ri_test_results$ri_lrt_p,
    ri_lrt_stat = ri_test_results$ri_lrt_stat,
    ri_var_subject = ri_test_results$ri_var_subject,
    ri_var_residual = ri_test_results$ri_var_residual,
    ri_icc = ri_test_results$ri_icc,
    # Original results
    n_obs = nrow(analysis_data),
    n_patients = length(unique(analysis_data$PatientID)),
    treatment_estimate = treatment_effect[1],
    treatment_se = treatment_effect[2],
    treatment_t = treatment_t,  # ä¿®æ­£: æ­£ã—ã„tå€¤ã‚’ä½¿ç”¨
    treatment_p = treatment_p,
    effect_size = effect_size,  # æ–¹å‘èª¿æ•´æ¸ˆã¿ã®Cohen's d
    effect_size_ci_lower = ci_lower_d,  # Cohen's d 95% CI lower
    effect_size_ci_upper = ci_upper_d,  # Cohen's d 95% CI upper
    icc = icc,
    residual_se = residual_sd,
    ci_lower = ci_result$lower,
    ci_upper = ci_result$upper,
    model_object = primary_model,
    analysis_data = analysis_data
  )
  
  # Add time delta results
  result_list <- c(result_list, time_delta_results)
  
  return(result_list)
}

# Function to run crossover analysis (fallback if mixed effects fail)
run_crossover_analysis <- function(data, metric) {
  metric_data <- data %>% filter(Evaluation == metric)
  
  cat(sprintf("\n=== Crossover Analysis for %s ===\n", metric))
  cat(sprintf("N observations: %d\n", nrow(metric_data)))
  
  complete_data <- metric_data %>% 
    filter(!is.na(Score_transformed) & !is.na(Baseline_z))
  
  if(nrow(complete_data) < 10) {
    cat("Insufficient data for analysis\n")
    return(NULL)
  }
  
  # Create difference scores
  patient_diffs <- complete_data %>%
    group_by(PatientID) %>%
    filter(n() == 2) %>%  # Only complete crossover data
    summarise(
      Score_diff = Score_transformed[Treatment == "aDBS"] - Score_transformed[Treatment == "cDBS"],
      Baseline_z_mean = mean(Baseline_z),
      .groups = 'drop'
    )
  
  if(nrow(patient_diffs) < 4) {
    cat("Insufficient crossover data\n")
    return(NULL)
  }
  
  # Simple t-test on difference scores
  t_test_result <- t.test(patient_diffs$Score_diff)
  
  # Model with baseline adjustment
  model_baseline <- lm(Score_diff ~ Baseline_z_mean, data = patient_diffs)
  model_summary <- summary(model_baseline)
  
  # Treatment effect is the intercept
  treatment_effect <- model_summary$coefficients["(Intercept)", ]
  
  # Calculate Cohen's D and apply direction adjustment
  effect_size_raw <- treatment_effect[1] / sigma(model_baseline)
  effect_size <- adjust_cohens_d_for_direction(metric, effect_size_raw)
  
  # Calculate 95% CI for Cohen's D (using accurate standard error formula)
  # Baseline adjusted crossover
  cohens_d_df <- df.residual(model_baseline)
  
  # Accurate SE for Cohen's d in paired/crossover design
  # SE(d) = sqrt(1/n + d²/(2*df)) @Hedges & Olkin (1985)
  n_patients <- nrow(patient_diffs)
  se_cohens_d <- sqrt(1/n_patients + (effect_size_raw^2)/(2*cohens_d_df))
  t_critical <- qt(0.975, cohens_d_df)
  
  # Raw Cohen's d CI (before direction adjustment)
  ci_lower_d_raw <- effect_size_raw - t_critical * se_cohens_d
  ci_upper_d_raw <- effect_size_raw + t_critical * se_cohens_d
  
  # Cohen's d CI after direction adjustment
  ci_lower_d <- adjust_cohens_d_for_direction(metric, ci_lower_d_raw)
  ci_upper_d <- adjust_cohens_d_for_direction(metric, ci_upper_d_raw)
  
  # Check CI order
  if(ci_lower_d > ci_upper_d) {
    temp <- ci_lower_d
    ci_lower_d <- ci_upper_d
    ci_upper_d <- temp
  }
  
  # For crossover analysis, we can't easily calculate time differences
  # because we don't have separate group estimates
  time_delta_results <- list(
    time_difference = NA,
    time_diff_ci_lower = NA,
    time_diff_ci_upper = NA,
    aDBS_time_adjusted = NA,
    cDBS_time_adjusted = NA,
    intuitive_ratio = NA,
    percent_reduction = NA,
    ratio_ci_lower = NA,
    ratio_ci_upper = NA
  )
  
  cat("\nTreatment Effect Summary (Crossover):\n")
  cat(sprintf("aDBS vs cDBS estimate: %.4f\n", treatment_effect[1]))
  cat(sprintf("p-value: %.4f\n", treatment_effect[4]))
  cat(sprintf("Effect size (Cohen's d, aDBS-favorable direction): %.4f [95%% CI: %.4f, %.4f]\n", effect_size, ci_lower_d, ci_upper_d))
  
  result_list <- list(
    metric = metric,
    analysis_type = "Crossover",
    carryover_p = NA,
    # Random Intercept LRT results (NA for crossover analysis)
    ri_lrt_p = NA,
    ri_lrt_stat = NA,
    ri_var_subject = NA,
    ri_var_residual = NA,
    ri_icc = NA,
    # Original results
    n_obs = nrow(patient_diffs) * 2,
    n_patients = nrow(patient_diffs),
    treatment_estimate = treatment_effect[1],
    treatment_se = treatment_effect[2],
    treatment_t = treatment_effect[3],
    treatment_p = treatment_effect[4],
    effect_size = effect_size,  # æ–¹å‘èª¿æ•´æ¸ˆã¿ã®Cohen's d
    effect_size_ci_lower = ci_lower_d,
    effect_size_ci_upper = ci_upper_d,
    icc = NA,
    residual_se = sigma(model_baseline),
    ci_lower = t_test_result$conf.int[1],
    ci_upper = t_test_result$conf.int[2],
    model_object = model_baseline
  )
  
  # Add time delta results
  result_list <- c(result_list, time_delta_results)
  
  return(result_list)
}

# Test if mixed effects models work
test_data <- ancova_data %>% 
  filter(Evaluation == "ADL") %>%
  filter(!is.na(Score_transformed) & !is.na(Baseline_z))

mixed_effects_available <- tryCatch({
  test_model <- lmer(Score_transformed ~ Treatment + (1|PatientID), 
                     data = test_data, REML = TRUE)
  TRUE
}, error = function(e) {
  cat("Mixed effects models failed:", e$message, "\n")
  FALSE
})

# Choose analysis method
if(mixed_effects_available) {
  cat("=== Using Mixed Effects Models ===\n")
  analysis_function <- run_mixed_effects_analysis
} else {
  cat("=== Using Crossover Analysis (Mixed Effects Failed) ===\n")
  analysis_function <- run_crossover_analysis
}

# Calculate Cohen's D and apply direction adjustment
cat("\n=== Cohen's D Interpretation Note ===\n")
cat("Cohen's d values are adjusted so that positive values always indicate aDBS superiority:\n")
cat("- Lower-the-better metrics (UPDRS, ADL, PDQ-39, dyskinesia/OFF time): negative treatment effect -> positive Cohen's d\n")
cat("- Higher-the-better metrics (ON time, MMSE): positive treatment effect â†’ positive Cohen's d\n\n")

# Run analysis for all metrics
ancova_results <- map(eval_items_ordered, ~analysis_function(ancova_data, .x)) %>%
  set_names(eval_items_ordered) %>%
  discard(is.null)

# Compile results into a summary table (UPDATED with Random Intercept LRT AND time delta columns)
results_summary <- map_dfr(ancova_results, function(result) {
  if(is.null(result)) return(NULL)
  
  tibble(
    Evaluation = result$metric,
    Analysis_Type = result$analysis_type,
    N_Observations = result$n_obs,
    N_Patients = result$n_patients,
    Carryover_P = result$carryover_p,
    # Random Intercept LRT columns (NEW)
    RI_LRT_P = ifelse(is.null(result$ri_lrt_p), NA_real_, result$ri_lrt_p),
    RI_LRT_Stat = ifelse(is.null(result$ri_lrt_stat), NA_real_, result$ri_lrt_stat),
    RI_Var_Subject = ifelse(is.null(result$ri_var_subject), NA_real_, result$ri_var_subject),
    RI_Var_Residual = ifelse(is.null(result$ri_var_residual), NA_real_, result$ri_var_residual),
    RI_ICC = ifelse(is.null(result$ri_icc), NA_real_, result$ri_icc),
    # Original columns
    Treatment_Estimate = result$treatment_estimate,
    Treatment_SE = result$treatment_se,
    Treatment_T = result$treatment_t,
    Treatment_P_uncorrected = result$treatment_p,
    CI_Lower = result$ci_lower,
    CI_Upper = result$ci_upper,
    Effect_Size = result$effect_size,
    Effect_Size_CI_Lower = ifelse(is.null(result$effect_size_ci_lower), NA_real_, result$effect_size_ci_lower),
    Effect_Size_CI_Upper = ifelse(is.null(result$effect_size_ci_upper), NA_real_, result$effect_size_ci_upper),
    ICC = result$icc,
    # UPDATED TIME COLUMNS with Delta Method
    Time_Difference = result$time_difference,
    Time_Diff_CI_Lower = result$time_diff_ci_lower,
    Time_Diff_CI_Upper = result$time_diff_ci_upper,
    aDBS_Time_Adjusted = result$aDBS_time_adjusted,
    cDBS_Time_Adjusted = result$cDBS_time_adjusted,
    Intuitive_Ratio = result$intuitive_ratio,
    Percent_Reduction = result$percent_reduction,
    Ratio_CI_Lower = result$ratio_ci_lower,
    Ratio_CI_Upper = result$ratio_ci_upper
  )
})

# FDR CORRECTION FOR TREATMENT EFFECTS ONLY

if(nrow(results_summary) > 0) {
  cat("\n=== APPLYING FDR CORRECTION (Benjamini-Hochberg) ===\n")
  
  # Extract uncorrected p-values for treatment effects
  treatment_p_values <- results_summary$Treatment_P_uncorrected
  
  # Apply Benjamini-Hochberg FDR correction
  treatment_p_fdr <- p.adjust(treatment_p_values, method = "BH")
  
  # Add FDR-corrected p-values to results
  results_summary_fdr <- results_summary %>%
    mutate(
      Treatment_P_FDR = treatment_p_fdr,
      Significant_uncorrected = Treatment_P_uncorrected < 0.05,
      Significant_FDR = Treatment_P_FDR < 0.05
    ) %>%
    # Reorder columns for clarity (including new RI LRT columns and time delta columns)
    select(Evaluation, Analysis_Type, N_Observations, N_Patients, 
           Carryover_P, 
           # Random Intercept LRT columns
           RI_LRT_P, RI_LRT_Stat, RI_Var_Subject, RI_Var_Residual, RI_ICC,
           # Treatment effect columns (CI placed immediately after Estimate for clarity)
           Treatment_Estimate, CI_Lower, CI_Upper, Treatment_SE, Treatment_T,
           Treatment_P_uncorrected, Treatment_P_FDR, 
           Significant_uncorrected, Significant_FDR,
           Effect_Size, Effect_Size_CI_Lower, Effect_Size_CI_Upper, ICC,
           # Time delta columns
           Time_Difference, Time_Diff_CI_Lower, Time_Diff_CI_Upper,
           aDBS_Time_Adjusted, cDBS_Time_Adjusted, Intuitive_Ratio, 
           Percent_Reduction, Ratio_CI_Lower, Ratio_CI_Upper)
  
  cat(sprintf("Number of tests corrected: %d\n", length(treatment_p_values)))
  cat("FDR correction applied using Benjamini-Hochberg method\n")
  
  print("\n=== STANDARD CROSSOVER ANALYSIS RESULTS (WITH FDR CORRECTION, RANDOM INTERCEPT LRT, AND DELTA METHOD) ===")
  print(results_summary_fdr)
  
  # Random Intercept LRT Summary (NEW)
  cat("\n=== RANDOM INTERCEPT LRT SUMMARY ===\n")
  ri_summary <- results_summary_fdr %>%
    filter(!is.na(RI_LRT_P)) %>%
    mutate(
      RI_Significant = RI_LRT_P < 0.05,
      RI_Status = case_when(
        RI_LRT_P < 0.01 ~ "Highly Significant",
        RI_LRT_P < 0.05 ~ "Significant", 
        RI_LRT_P < 0.10 ~ "Marginal",
        TRUE ~ "Non-significant"
      )
    ) %>%
    select(Evaluation, RI_LRT_P, RI_LRT_Stat, RI_Var_Subject, RI_Var_Residual, RI_ICC, RI_Status) %>%
    arrange(RI_LRT_P)
  
  print(ri_summary)
  
  cat(sprintf("\nRandom intercept significant (p < 0.05): %d/%d metrics\n", 
              sum(ri_summary$RI_LRT_P < 0.05, na.rm = TRUE), 
              nrow(ri_summary)))
  
  time_variables_results <- results_summary_fdr %>%
    filter(str_detect(Evaluation, "time")) %>%
    filter(!is.na(Time_Difference))
  
  if(nrow(time_variables_results) > 0) {
    cat("\n=== TIME VARIABLES: DELTA METHOD RESULTS ===\n")
    print(time_variables_results %>%
            select(Evaluation, Treatment_P_FDR, Significant_FDR, 
                   Time_Difference, Time_Diff_CI_Lower, Time_Diff_CI_Upper,
                   aDBS_Time_Adjusted, cDBS_Time_Adjusted, 
                   Intuitive_Ratio, Percent_Reduction))
    
    cat("\n=== TIME VARIABLES: CLINICAL INTERPRETATION WITH DELTA METHOD ===\n")
    time_variables_results %>%
      mutate(
        Clinical_Interpretation = case_when(
          str_detect(Evaluation, "OFF|Dys") & Time_Difference < 0 ~ 
            paste0("aDBS reduced ", Evaluation, " by ", round(abs(Time_Difference), 2), " hours (", 
                   round(abs(Percent_Reduction), 1), "% improvement)"),
          str_detect(Evaluation, "OFF|Dys") & Time_Difference > 0 ~ 
            paste0("aDBS increased ", Evaluation, " by ", round(Time_Difference, 2), " hours (", 
                   round(Percent_Reduction, 1), "% worsening)"),
          str_detect(Evaluation, "ON") & Time_Difference > 0 ~ 
            paste0("aDBS increased ", Evaluation, " by ", round(Time_Difference, 2), " hours (", 
                   round(abs(Percent_Reduction), 1), "% improvement)"),
          str_detect(Evaluation, "ON") & Time_Difference < 0 ~ 
            paste0("aDBS reduced ", Evaluation, " by ", round(abs(Time_Difference), 2), " hours (", 
                   round(Percent_Reduction, 1), "% worsening)"),
          TRUE ~ "No clear interpretation"
        ),
        Delta_CI_Interpretation = paste0("95% CI: [", round(Time_Diff_CI_Lower, 2), 
                                         ", ", round(Time_Diff_CI_Upper, 2), "] hours")
      ) %>%
      select(Evaluation, Clinical_Interpretation, Delta_CI_Interpretation) %>%
      {
        for(i in seq_len(nrow(.))) {
          cat(paste0("  • ", .$Evaluation[i], ": ", .$Clinical_Interpretation[i], 
                     " (", .$Delta_CI_Interpretation[i], ")\n"))
        }
      }
  }
  
  # Analysis type summary (all metrics use full crossover analysis)
  cat("\n=== ANALYSIS TYPE SUMMARY ===\n")
  cat(sprintf("All metrics analyzed using full crossover analysis: %d metrics\n", nrow(results_summary_fdr)))
  
  # Significant results AFTER FDR correction
  significant_results_fdr <- results_summary_fdr %>%
    filter(Significant_FDR == TRUE) %>%
    arrange(Treatment_P_FDR)
  
  cat("\n=== Significant Treatment Effects AFTER FDR Correction (q < 0.05) ===\n")
  if(nrow(significant_results_fdr) > 0) {
    print(significant_results_fdr %>% 
            select(Evaluation, Analysis_Type, Treatment_Estimate, 
                   Treatment_P_uncorrected, Treatment_P_FDR, 
                   CI_Lower, CI_Upper, Effect_Size))
  } else {
    cat("No significant treatment effects after FDR correction.\n")
  }
  
  # Results that were significant before but not after FDR correction
  lost_significance <- results_summary_fdr %>%
    filter(Significant_uncorrected == TRUE & Significant_FDR == FALSE) %>%
    arrange(Treatment_P_uncorrected)
  
  cat("\n=== Results That Lost Significance After FDR Correction ===\n")
  if(nrow(lost_significance) > 0) {
    print(lost_significance %>% 
            select(Evaluation, Analysis_Type, Treatment_Estimate,
                   Treatment_P_uncorrected, Treatment_P_FDR, Effect_Size))
  } else {
    cat("No results lost significance after FDR correction.\n")
  }
  
  # Marginally significant results after FDR correction
  marginal_results_fdr <- results_summary_fdr %>%
    filter(Treatment_P_FDR >= 0.05 & Treatment_P_FDR < 0.10) %>%
    arrange(Treatment_P_FDR)
  
  cat("\n=== Marginally Significant Treatment Effects After FDR (0.05 <= q < 0.10) ===\n")
  if(nrow(marginal_results_fdr) > 0) {
    print(marginal_results_fdr %>% 
            select(Evaluation, Analysis_Type, Treatment_Estimate,
                   Treatment_P_uncorrected, Treatment_P_FDR, 
                   CI_Lower, CI_Upper, Effect_Size))
  } else {
    cat("No marginally significant treatment effects after FDR correction.\n")
  }
  
  cat("\n=== Effect Size Summary (Positive = aDBS Favorable) ===\n")
  effect_size_summary <- results_summary_fdr %>%
    arrange(desc(Effect_Size)) %>%
    select(Evaluation, Effect_Size, Treatment_P_FDR, Significant_FDR)
  
  print(effect_size_summary)
  
  # Summary of carryover effects (unchanged - no correction applied)
  cat("\n=== DIFFERENTIAL CARRYOVER SUMMARY (no correction applied) ===\n")
  carryover_summary <- results_summary_fdr %>%
    filter(!is.na(Carryover_P)) %>%
    mutate(
      Carryover_Status = case_when(
        Carryover_P < 0.05 ~ "Significant",
        Carryover_P < 0.10 ~ "Marginal", 
        TRUE ~ "Non-significant"
      )
    ) %>%
    select(Evaluation, Carryover_P, Carryover_Status, Analysis_Type)
  
  print(carryover_summary)
  
  # Delta Method Summary for Time Variables
  if(nrow(time_variables_results) > 0) {
    cat("\n=== DELTA METHOD SUMMARY FOR TIME VARIABLES ===\n")
    delta_summary <- time_variables_results %>%
      mutate(
        CI_Width = Time_Diff_CI_Upper - Time_Diff_CI_Lower,
        CI_Contains_Zero = (Time_Diff_CI_Lower <= 0) & (Time_Diff_CI_Upper >= 0),
        Precision = case_when(
          CI_Width < 1 ~ "High precision (CI < 1 hour)",
          CI_Width < 2 ~ "Moderate precision (CI < 2 hours)",
          TRUE ~ "Low precision (CI >= 2 hours)"
        )
      ) %>%
      select(Evaluation, Time_Difference, CI_Width, CI_Contains_Zero, Precision)
    
    print(delta_summary)
    
    cat("\nDelta Method Notes:\n")
    cat("- CI_Contains_Zero = TRUE indicates time difference not significantly different from zero\n")
    cat("- Narrower CI widths indicate more precise estimates\n")
    cat("- All calculations use Delta method for non-linear transformations\n")
  }
  
  # Export results with FDR correction, Random Intercept LRT, AND DELTA METHOD
  write_csv(results_summary_fdr, "results_standard_crossover.csv")
  write_csv(carryover_summary, "results_carryover_effects.csv")
  write_csv(ri_summary, "results_random_intercept_lrt.csv")
  write_csv(ancova_data, "data_ancova_long_format.csv")
  saveRDS(ancova_results, "results_standard_crossover.rds")
  
  # Export time variables summary separately with Delta Method results
  if(nrow(time_variables_results) > 0) {
    write_csv(time_variables_results, "results_time_variables_delta_method.csv")
    write_csv(delta_summary, "results_delta_method_summary.csv")
  }
  
  cat("\n=== Files Exported ===\n")
  cat("1. Standard_Crossover_Results_Complete_Final.csv - Complete analysis results with ALL features\n")
  cat("2. Carryover_Effects_Summary.csv - Differential carryover summary\n")
  cat("3. Random_Intercept_LRT_Summary.csv - Random intercept LRT summary\n")
  cat("4. ANCOVA_long_format_data.csv - Transformed data\n")
  cat("5. Standard_Crossover_Analysis_Complete_Final.rds - Complete analysis objects\n")
  if(nrow(time_variables_results) > 0) {
    cat("6. Time_Variables_DeltaMethod_Results.csv - Time variables with Delta method CIs\n")
    cat("7. Delta_Method_Summary.csv - Delta method precision summary\n")
  }
  
} else {
  cat("No successful analyses completed.\n")
}

cat("\n=== FINAL ANALYSIS COMPLETE ===\n")
cat("Complete crossover analysis implemented with ALL advanced features:\n")
cat("1. Differential carryover testing (Treatment x Period interaction via LRT)\n")
cat("2. Random intercept LRT testing (subject variance significance)\n")
cat("3. Adaptive analysis strategy (full crossover vs first period only)\n")
cat("4. FDR correction (Benjamini-Hochberg) for multiple comparisons\n")
cat("5. Cohen's d with direction adjustment (positive = aDBS favorable)\n")
cat("6. Delta method for rigorous time difference estimation\n")
cat("7. Comprehensive variance component analysis (ICC, subject/residual variance)\n")
cat("8. Clinical interpretation with absolute time differences\n")
cat("9. Precision assessment for all estimates\n")
cat("10. Complete export of all results and summaries\n")